clear;
n = 64;
W = [1, 2/3];
figure('units', 'normalized', 'position', [0, 0, 0.5, 1]);
for j=1:2
    clf;
    w = W(j);
   
    T = zeros(n-1);
    for i=1:n-1
        if i ~= n-1
            T(i, i+1) = w/2;
            T(i+1, i) = w/2;
        end
        T(i, i) = 1 - w;
    end

    for k=1:n-1
        num=0;
        is_stop = false;
        deter = true;
        for i = 1:n-1
            X(i, 1) = sin(i/n*k*pi);
        end
        while ~is_stop
            num = num + 1;
            X=T*X;
            deter = true;
            for i = 1:n-1
                if abs(X(i)) > 1/100*abs(sin(1/n*k*pi))
                    deter = false;
                end
            end
            if deter
                is_stop = true;
            end
            if num == 100
                is_stop = true;
            end
        end
        number(k) = num;
    end
    K = 1:1:n-1;
    plot(K, number, '-k');
    hold on;
    axis([-10, n, 0, 110]);
    plot([0, n], [0, 0], '-k');
    plot([0, 0], [0, 100], '-k');
    plot([16, 16], [0, 50], '--k');
    if j== 1
        plot([48, 48], [0, 50], '--k');
    end
    text(n/2, -10, 'Wavenumber k', 'Fontsize', 25, 'color', 'k');
    text(-5, 50, 'Iterations', 'Fontsize', 25, 'color', 'k', 'Rotation', 90);
    if j == 1
        s = 'w=1';
    else
        s = 'w=2/3';
    end
    text(n/2, 50, s, 'Fontsize', 25, 'color', 'k');
    for i = 1:8
        plot([n*i/8, n*i/8], [0, 1], '-k');
        plot([0, 0.5], [100*i/8, 100*i/8], '-k');
        if mod(i, 2) == 0
            text(n*i/8-1, -1.5, num2str(n*i/8), 'Fontsize', 25, 'color', 'k');
            l=0;
            if i == 8
                l = -0.5;
            end
            text(-2.5 + l, 100*i/8, num2str(100*i/8), 'Fontsize', 25, 'color', 'k');
        end
    end
    plot([0, n], [number(16), number(16)], '--k');
    axis off;
    s = ['../Exercise9_28_', num2str(j), '.eps'];
    print(gcf,'-depsc', s);
    hold off;
end